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Abstract 

Motivated by recent experimental observations, we study theoretically multiple bright solitary waves of trapped Bose- 
Einstein condensates. Through variational and numerical analyses, we determine the threshold for collapse of these 
states. Under 7r-phase differences between adjacent waves, we show that the experimental states lie consistently at the 
threshold for collapse, where the corresponding in-phase states are highly unstable. Following the observation of two 
long-lived solitary waves in a trap, we perform detailed three-dimensional simulations which confirm that in-phase 
waves undergo collapse while a 7r-phase difference preserves the long-lived dynamics and gives excellent quantitative 
agreement with experiment. Furthermore, intermediate phase differences lead to the growth of population asymmetries 
between the waves, which ultimately triggers collapse. 
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1. Introduction 

The nonlinearity present in atomic Bose-Einstein 
condensates has led to demonstrations of fascinat- 
ing self-trapped states known as solitons. These one- 
dimensional wavepackets, well known in non-linear 
optics and other fields pQ, arise when the nonlin- 
earity of the medium counter-acts the effects of dis- 
persion. They have been realized in several distinct 
matter- wave forms: bright |2I3I4| . dark [5] and gap 
6 solitons. The bright solitons manifest themselves 
as self-trapped lumps of matter, held together by 
attractive atomic interactions. In three-dimensions, 
bright matter-wave "solitons" are self-trapped in 
one dimension and require external confinement in 
the remaining two directions. We will henceforth 
refer to these 3D solitonic states as bright solitary 
waves (BSWs). Due to their self-trapping proper- 
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ties, matter-wave BSWs offer significant possibili- 
ties in atom-optical applications such as atom in- 
terferometry [2j and probing the atom-surface in- 
teraction [7]. However, the three-dimensional na- 
ture of BSWs leads to the presence of an undesir- 
able collapse instability when the attractive interac- 
tions become too strong [8 9 10 11 . Not only does 
this affect the static properties of BSWs, it can also 
destabilise their collisions |12I13I14|T5] . Although 
techniques to suppress collapse effects in attractive 
BECs have been proposed, e.g., by applying rotation 
about the z-axis [16] or by time-modulating the scat- 
tering length [T7] , the collapse instability ultimately 
remains. A detailed understanding of the properties 
of matter- wave BSWs and their regimes of collapse is 
therefore essential for the advancement of this field. 

The presence of solitonic solutions is revealed by 
considering the mean-field and zero-temperature 
limit of atomic BECs. Here the BEC 'wavefunction' 
ip{v, t) satisfies a nonlinear wave equation known as 



Preprint submitted to Elsevier 



28 February 2008 



the Gross-Pitaevskii equation (GPE) [TB] , given by, 

^ (1) 



where m is the atomic mass. The nonlinear term 
arises from the short-range atomic interactions, 
characterised by the s-wave scattering length a s , 
and can be repulsive (et s > 0) or attractive (a s < 0). 
In ID and in the absence of an external potential 
Kxt(r) = 0, this has an identical form to the ID cu- 
bic nonlinear Schrodinger equation and supports the 
exact dark and bright soliton solutions derived by 
Zakharov and Shabat [19] . Key properties of these 
ID solitons are that they can have any population, 
their collisions are elastic [5U] and they are stable 
to thermal dissipation. However, atomic BECs are 
intrinsically 3D objects, typically confined by har- 
monic traps of the form V ext (r) = moj 2 (r 2 +X 2 z 2 )/2, 
where uj r and \uj r are the radial and axial trap 
frequencies, respectively. The presence of the extra 
dimensions modifies the special properties of the ID 
soliton by introducing a critical atomic population, 
inelastic collisions [15] and thermal dissipation [21] . 

The critical population and inelastic collisions of 
BSWs arise due to the collapse instability which af- 
fects attractively-interacting BECs in general. It is 
convenient to introduce the dimensionless interac- 
tion parameter k defined as, 



k = 



N\a s 



(2) 



where N is the number of atoms in the system and 
a r = y^h/mujr is the radial harmonic oscillator 
length of the trap0. Collapse of the system occurs 
when k exceeds a critical value k c , which is typi- 
cally of the order of unity. This implies that there 
is a critical population N c beyond which collapse 
occurs. Note that it is the presence of trapping that 
leads to a finite value of N c , whereas a homogeneous, 
untrapped BEC is always unstable to collapse ,26] . 

In a recent experiment multiple BSWs were gen- 
erated in a three-dimensional trap [4] . Intriguingly, 
the ensuing dynamics were remarkably robust and 
long-lived. In this work we analyse the properties 
of bright solitary matter-waves under external con- 
fining potentials with particular emphasis on this 
experiment [4 1. After reviewing the details of this 
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experiment (Section 2), we employ a variational 
approach and full solution of the Gross-Pitaevskii 
equation to study the ground state and first-excited 
state solutions of the system (Section 3). We extend 
this to additional solitary waves using a dynamical 
model, revealing the threshold of the collapse in- 
stability for up to four solitary waves and compare 
to experimental measurements. We then directly 
simulate the experimental oscillations of two BSWs 
(Section 4) and show that a 7r-phase difference is es- 
sential to support the observed dynamics and gives 
excellent agreement with the experimental data. 
Finally, we present the conclusions of our work 
(Section 5). 



2. The JILA experiment 

We will review the BSW experiment of Cornish 
et al. 4 and its observations. Firstly, a stable 85 Rb 
BEC was formed with repulsive s-wave interactions. 
This typically contained 15000 atoms with less than 
500 thermal atoms. The magnetic confining trap was 
cylindrically symmetric with radial frequency Lu r = 
2-7T x 17.3Hz and is weakly elongated with a trap ra- 
tio A = 0.4. For this system, experimental [25] and 
theoretical work [ll|23j agrees that the critical inter- 
action parameter for collapse is k c 0.64. Follow- 
ing previous experiments [25] , the s-wave scattering 
length was then quickly tuned to be attractive by 
means of a molecular Feshbach resonance [37] ■ At 
this point the number of atoms in the system greatly 
exceeded the critical number triggering a collapse. 
During the collapse, three-body atomic losses rise 
and eventually stabilised the condensate. The "rem- 
nant" condensate typically contained more than the 
critical number of atoms, N c , and was clearly di- 
vided in the axial direction into a symmetric ar- 
rangement of distinct wavepackets, i.e. bright soli- 
tary waves, that oscillated along the weaker axial 
direction of confining potential maintaining a sta- 
tionary centre of mass. Note that it is thought that 
the collapsing condensate fragments into multiple 
wavepackets via a modulational instability [12128] , 
Up to six BSWs were observed, depending on fac- 
tors such as the final scattering length and the ini- 
tial number of atoms [4] . In Fig. Q] we present the 
typical appearance of the condensate column den- 
sity following the formation of the BSWs. The im- 
ages were taken when the BSWs reached the outer 
turning points of their oscillatory motion in the har- 
monic trap and show a cross-section of the optical 
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Fig. 1. Axial profile of the optical depth from the JILA ex- 
periment following BSW formation showing (a) two BSWs 
and (b) three BSWs. The profiles were taken when the BSWs 
reached the outer turning points of their oscillatory motion 
in the harmonic trap. In each case, the profile is fitted with 

the form ^=1™ [ a i exp {-(z - ft) 2 /^ 2 }] , where N BSW 
is the number of BSWs and ai , ft and ai arc fitting param- 
eters. In (a), the mean BSW width is "a ~ 7.8(4)/xm and in 
(b) it is a a 4.5(4)/im. 

depth in the axial direction for the case of two and 
three BSWs. Note that the density profiles are ap- 
proximately symmetric about the origin. In particu- 
lar, quantitative measurements of the experimental 
system were made for fixed scattering length a s = 
— 0.6nm and approximately 4000 atoms. Since the 
system initially contained approximately 500 ther- 
mal atoms, we will assume the number of conden- 
sate atoms to be TV = 3500. This corresponds to 
an interaction parameter of k — 0.8 which exceeds 
the critical interaction parameter k c . Despite this, 
the observed dynamics were surprisingly stable with 
negligible dissipation over 3 s. For the case of two 
BSWs, the dynamics consist of oscillations in anti- 
phase along the axial direction of the trap, with the 
BSWs colliding repeatedly at the trap centre. It is 
thought that the stable dynamics were supported by 
the existence of a repulsive 7r-phase difference be- 
tween the BSWs, which allowed each BSW to sup- 
port an atom number corresponding to an interac- 
tion parameter just below k c . Indeed, experimental 
measurements for up to four BSWs showed that the 
the average value of k per BSW never exceeded k c 
and this has been recently found to be in good agree- 
ment with a theoretical study [TT]. Note that the 
observed BSW dynamics showed no significant ther- 
mal dissipation, despite the presence of the highly 
energetic burst of atoms ejected from the conden- 
sate during the collapse [59]. It is likely that these 
excitations are so "hot" (Tburst ~ 50 nK) and dilute 
that thermal equilibrium is not reached over the ex- 
perimental timescales, with the "thermal atoms" re- 
maining effectively invisible to the BSWs through- 
out. Note that the existence of a 7r-phase difference 



has also been inferred in the experiment of Strecker 
et al. 3^ by the observed repulsive interactions be- 
tween the BSWs. 

3. BSW solutions and the collapse instability 

In 3D and in the presence of interactions there are 
no exact analytic solutions of the GPE and solutions 
must be obtained numerically or via approximated 
approaches. In the latter case, variational methods 
[8I9I10I1I] provide considerable insight without the 
the need for full numerical solutions. Here we will 
consider both a variational approach and numerical 
solution of the GPE. 

3.1. Ground state solutions 

Firstly, we shall study the ground state solutions 
in the system. We approximate the ground state 
solutions (denoted by the 0-subscript) by a single- 
peaked ansatz of the form [819110111] , 



N 



where L z and L r represent the axial and radial sizes, 
respectively. The energy of the system e is defined 
by the GP energy functional via, 



e = d 6 r 



r ^2 



2m 



IW>| 2 + K xt H : 



M 4 



■(4) 



By substituting the ansatz into the energy func- 
tional of Eq. ([4]) we determine the variational energy 
for the ground state Eo- For convenience we employ 
the rescaled parameters l r = L r /a r , l z = L z /a r and 
E = e/(Nhw r ), where a r = \J h) muj r is the radial 
harmonic oscillator length. Furthermore, we intro- 
duce our interaction parameter k — N\a s \/a r . The 
ground state ansatz energy then becomes, 



1/1 1 



i" 2 



X 2 li 



k 



2nl z l 2 



.(5) 



The first group of terms represents kinetic energy, 
the second group represents potential energy and 
the final term represents the energy arising from the 
s-wave interatomic interactions. This equation de- 
fines an energy landscape for the system in terms of 
l z and l r , and variational solutions correspond to lo- 
cal energy minima in this landscape. We denote the 
widths of the variational solution by l z and 1®. The 
variational method has been employed successfully 
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to study the ground states of the system [8111112] . 
In particular, in the absence of axial trapping (A = 
0), local energy minima arise which correspond to 
self-trapped BSW solutions. When k exceeds a crit- 
ical value the local energy minimum ceases to exist 
and the global energy minimum, which occurs at the 
origin, dominates the system. This corresponds to 
collapse. 

We have calculated the variational solutions un- 
der an axial trap defined by A = 0.4. The size 
of the variational solutions as a function of in- 
teraction parameter k are presented in Fig. HJa). 
For k = 0, the solution corresponds to the exact 
non-interacting gaussian ground state with l z = 
y/h/mujz — A -1 / 2 a r = 1.58a r . As the attractive 
interactions grow in size, the solutions, which are 
always elongated in the z— direction, shrink in both 
dimensions. Finally, at k = 0.746, the variational 
solutions disappear and the system is unstable to 
collapse. 

We have also calculated the exact solutions by 
solving the GPE numerically. This is performed us- 
ing the imaginary-time propagation technique: un- 
der the substitution t — > —it in the GPE, the equa- 
tion evolves to the ground state of the system, pro- 
viding it exist£El The axial and radial lengthscales 
of the GPE solutions correspond to the distance over 
which the density decreases by a factor 1/e from 
its peak value. The GPE lengthscales are presented 
in Fig. [UJa) by crosses (filled circles) for the radial 
(axial) direction. For k = the GPE results agree 
exactly with the variational prediction. As k is in- 
creased the GPE predictions decrease at a faster 
rate than the variational method. Furthermore, the 
GPE solutions become unstable to collapse for k c = 
0.630(5). This is approximately 15% lower than the 
variational method and arises because the varia- 
tional method consistently over-estimates the BSW 
widths and therefore under-estimates the peak den- 
sity, which is the trigger of collapse. This result is 
consistent with previous studies |ll|23j . where k c 
was mapped out for a range of trap ratios. In a pre- 
vious JILA experiment [25], the critical interaction 
parameter for collapse (of a single wavepacket) was 
measured to be k c = 0.64(7) [3 This is in excel- 
lent agreement with the GPE prediction of k c — 
0.630(5), as noted elsewhere |ll|23j . 



2 Note that the equation is no longer unitary and so must 
be renormalised at each time step to preserve atom number 

3 Note that the quoted value of k c in [25] was revised in |30| 
based on more accurate measurements of a B . 
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Fig. 2. (a) Ground state solutions under axial trapping 
A = 0.4 according to the variational approach of Eq. J5jl 
(solid/dashed lines) and by numerical solution of the GPE 
(crosses/circles). Upper and lower lines indicate the axial 
lengthscale 1® and radial lengthscale ij! of the solution, re- 
spectively, (b) First-excited state solutions for A = 0.4 ac- 
cording to the variational energy of Eq. J7J and the full nu- 
merical solution of the GPE. 

Comparison of the experimental BSW density 
profiles (e.g. Fig. [1]) to the theoretical predictions 
is, however, unsuccessful. According to the theo- 
retical results of Fig. [Ha), the axial width of the 
BSW should be at the very most equal to the non- 
interacting value of l° z (k = 0) = 1.6a r « 4//m. How- 
ever, the mean fitted width of each BSW is approx- 
imately 7.8/im in Fig.rjja) and 4.5/j,m in Fig.QJb). 
This discrepancy is almost certainly due to the low 
resolution of the experimental imaging [4j . 



3.2. First-excited state solutions 

Since two (and more) BSWs were also observed 
in the JILA experiment, and are believed to be sup- 
ported by a 7r-phase difference, it is pertinent to con- 
sider the first-excited state of the system (denoted 
by 1-subscript). We extend the variational approach 
by replacing the gaussian axial profile with that of 
the first-excited harmonic oscillator state, such that 
the ansatz is, 
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^ z)= i^hi {£) 

XeXP ("2Zf) eXP (-^)' (6) 

The density profile of this ansatz is double-peaked 
and can be interpreted as two BSWs featuring a tv- 
phase difference. Note that a similar approach by 
Michincl et al. 31j employed Hermite functions to 
study multiple BSWs in a trap. 

Following the same method as for the ground state 
ansatz, we arrive at the variational energy for the 
first-excited state, 

e- 1 ( 1 + 3 \ 1 / 2 ^xni \ 



We have obtained the corresponding first-excited 
variational solutions and plotted their lengthscales 
in Fig.[5{b). Collapse of the first-excited variational 
solution occurs at k c — 1.508. We have also ob- 
tained the exact first-excited states of the full GPE, 
shown in Fig. EJb) by crosses. This was performed 
using the imaginary-time technique while enforc- 
ing the wavefunction to be asymmetric via ip{ z ) — 
—ip(—z). These solutions become unstable at k c = 
1.145(5). Note that in the absence of axial confine- 
ment, there is no stationary first-excited state, be- 
cause the BSWs exert a repulsive force on each other 
which decays exponentially with their separation. 
The lowest energy state is therefore when the BSWs 
are infinitely separated. 

According to both the variational ansatz and the 
full GPE solution, the first-excited state supports an 
interaction parameter which is almost twice that of 
the ground state. Note that if the individual BSWs 
were completely independent we would expect the 
system to support exactly twice k c of the individual 
BSWs. 

3.3. Up to four BSWs 

In the JILA experiment, quantitative measure- 
ments were made of up to four BSWs, with the re- 
sults presented in Fig. [3] (points with error bars). 
While the total atom number typically exceeded 
the critical atom number for the ground state N c , 
the average atom number per BSW was less than, 
or approximately equal to, N c . This observation is 
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Fig. 3. The ratio Ntot/N c as a function of the num- 
ber of BSWs A^bsw- GPE simulations show the critical 
points for phases differences of A</> = (dashed line) and 
A<f> = 7r (bold solid line), for a fixed scattering length of 
a B = — 0.6nm. Above/below these lines the configurations 
are stable/unstable to collapse. The experimental JILA data, 
shown by points with error bars, is obtained at various scat- 
tering lengths. The corresponding GPE results are presented 
by black crosses, where it is assumed that 10% of the exper- 
imentally detected atoms were non-condensed. The function 
Atot = AbswA c is shown for comparison (grey line). 

thought to be a direct consequence of repulsive tt- 
phase differences between the BSWs. We now con- 
sider configurations of up to four BSWs. In prin- 
ciple one could extend the variational approach to 
model any number of BSWs by using higher-order 
excited states of the harmonic oscillator. However, 
here we will employ numerical solutions of the full 
GPE. States of one and two BSWs are obtained 
by imaginary time propagation as detailed above. 
States of three and four BSWs cannot be formed 
by imaginary time propagation since the lowest en- 
ergy state, which is where the atoms populate the 
central BSWs, is unstable to collapse. A dynami- 
cal method is employed |llj where we begin with a 
ground state repulsively-interacting BEG The in- 
teractions are then switched to the required attrac- 
tive value while simultaneously imposing a periodic 
distribution of 7r-phase steps. In this manner dy- 
namic states of three or four BSWs are created and 
their critical regime for collapse can be probed. The 
results for a fixed scattering length of a s = —0.6 nm 
are shown by the solid line in Fig. [3] Above (below) 
this line, the configurations are unstable (stable). 
We see that four BSWs are readily supported, pro- 
viding 7r-phasc differences are present. For 0-phase 
differences, only one BSW with up to A" c can be sup- 
ported, as indicated by the dashed line. Assuming 
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the BSWs to be independent and that each contains 
up to N c atoms, the function N to t — ^bsw-^c (grey 
line) is satisfied. The numerical results deviate from 
this function as iVasw increases. This is due to the 
presence of interactions between the BSWs and the 
unequal distribution of atoms in the BSWs, i.e. the 
central (outer) BSWs contain more (less) atoms, as 
observed in the experiment. 

The experimental number measurements (points 
with error bars) were obtained at various scattering 
lengths. Using these scattering lengths, we have nu- 
merically evaluated the critical points for collapse 
according to the GPE. Note that we have assumed 
10% of the experimental atom number to be non- 
condensed. The experimental data is in excellent 
agreement with these GPE predictions, with every 
data point showing consistency between theory and 
experiment. These results show the importance of 
7r-phasc differences. Furthermore, they show that 
the experimental system consistently forms a state 
which is right at the limit of collapse. This is a re- 
markable effect, given that the system is initially in 
a state which is highly unstable to collapse. 

4. Dynamics of two BSWs in a trap 

The collision of two BSWs in a homoge- 
neous waveguide has been considered previously 
[13114115] . During the collision, a high density state 
forms during the collapse and providing the in- 
teraction parameter is sufficiently large, this can 
induce collapse. The collisions are most prone to 
collapse when (i) the BSWs are in-phase, since this 
maximises the overlap and hence the peak density 
during the collision, and (ii) for low speed collisions, 
since the timescale over which this overlap occurs is 
large. In contrast, the collapse instability is heavily 
suppressed when the BSWs feature a 7r-phase dif- 
ference or if the impact speed is high, since this re- 
duces the timescale over which a collapse can occur. 
Furthermore, we have recently shown that when the 
relative phase A0 between the colliding BSWs does 
not equal zero or 7r, a sizeable population transfer 
can occur between the BSWs [15] , For < A0 < ir, 
this population transfer flows in one direction, while 
for 7r < A<p < 2tt, it flows in the opposite direction. 
While these studies have typically considered single 
collisions in a homogeneous waveguide, the presence 
of axial trapping enables multiple collisions. One- 
dimensional approaches to BSW collisions in a trap 
have been made, with a recent study employing a 



particle model for the BSWs [35). Below we make 
a detailed study of the 3D dynamics of two BSWs 
under the conditions of the JILA experiment. 

In the experiment, the axial density profile was 
measured at regular intervals and each profile was 
fitted to a single gaussian profile to give the ax- 
ial full- width- half-maximum (FWHM). The exper- 
imental data showing the evolution of the FWHM 
is presented in Fig.[4jiii). The key observations are 
that the FWHM oscillates primarily at 2\ui r and 
that there is negligible dissipation over 3 s which cor- 
responds to 40 oscillations. As discussed in Section^ 
the condensate contains approximately 3500 atoms. 
We assume each BSW to contain N = 1750 atoms 
and therefore k « 0.4. At the extreme points of their 
motion, the BSWs were observed to be displaced 
by approximately zq = 16/im from the trap centre. 
Note this this relatively large separation means that 
it is more appropriate here to consider the system 
as a collection of single BSWs rather than a true 
first-excited state of the system. Our initial state 
therefore consists of two ground state BSWs with 
k = 0.4, positioned at zq = ±16/im. In addition, 
we also impose a relative phase difference A(f> be- 
tween the initial BSWs. Typically, the experimental 
BSW density profile remained symmetric about the 
trap centre, as illustrated in the experimental pro- 
files shown in Fig. [T] Since only relative phase dif- 
ferences of and it (modulo 2tt) preserve this sym- 
metry in BSW collisions [15] we will concentrate on 
these phase differences. 

Note that the corresponding dynamics for three 
BSWs in a trap appear significantly less stable than 
the two-BSW state, even with 7r-phase differences. 
This will be considered in a future study. 

4.1. Relative phases A<p = and it 

Simulation of the dynamics for short times are 
shown in Fig. |4ja) f° r W ^ ( t ) — an d (ii) A<fr = 
ir. Initially the BSWs accelerate identically towards 
the trap centre. As they interact the relative phase 
becomes apparent with the BSWs overlapping for 
A(f> = and 'bouncing' for A(j) = n. However, apart 
from at the point of interaction, the BSW dynamics 
are almost identical over this time scale. The corre- 
sponding FWHM is presented in Fig.[4f a) (iii) and is 
in excellent agreement with the experimental data. 
The 2Xui r oscillations in the FWHM arises since each 
BSW oscillates in the axial direction at a frequency 
of Xu! r . Note that while the simulated FWHM be- 
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Fig. 4. (a) Evolution of the radially-integrated axial density 
at early times for two BSWs (JVg = 1750, zo = ±16^tm, 
uj z /2it = 6.8 Hz and u) r /2-K = 17.3 Hz) with (i) A<j> = 
and (ii) A(j> = n. (iii) Full-width-half-maximum FWHM of 
a gaussian fit to the axial density for A<j> = (dashed line), 
A<f> = tt (solid line), and experimental data (points) (4]. Note 
that in order to match the phase of the oscillations we have 
shifted the experimental data in time, (b) Same as (a) but 
for late times. 

gins at a maximum, the experimental data begins at 
a minimum [4] since the BSWs are created in close 
proximity. Consequently, we have shifted the exper- 
imental data in Fig. [4] by a quarter of the oscillation 
period such that the experimental and simulated 
data are in phase. Also note that due to the resolu- 
tion limit in the experiment, the FWHM cannot be 
resolved below approximately 15/im. 

The corresponding dynamics for late times are 
shown in Fig. H^b). The A<p = BSWs have fully 
collapsed, no longer matching the experimental re- 
sults. In contrast, the A(j) = tt collisions remains 
practically unaffected. Furthermore, they are in ex- 
cellent agreement with the experimental data. This 
shows that a phase difference of tt is crucial to sup- 
port the long-lived oscillations observed in the JILA 
experiment. The presence of 7r-phasc differences was 
also inferred in the BSW experiment of Strecker et 
al. [3]. 

4.2. Intermediate relative phases < A<j) < tt 

For comparison, we have also simulated the dy- 
namics of the two BSWs for intermediate relative 
phases in the range < A<fi < tt. The density dy- 
namics are presented in Fig. [5] for (a) A<fi — 7r/4, 
(b) 7r/2 and (c) 3tt/A. Additionally, in Fig.E^d) we 
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Fig. 5. Dynamics of two BSWs (k = 0.4) in the A = 0.4 JILA 
trap with (a) A<j> = tt/4, (b) A</> = tt/2 and (c) A</> = 3tt/4. 
The initial separation is 2q = 16^tm. (d) Population differ- 
ence across the system (TVl — Nji)/N for A<f> = 7r/4 (solid 
blue line), <j> = it/ 2 (dotted black line) and 3-7r/4 (dashed red 
line). For the cases of A<p = and it, the population differ- 
ence is zero at all times, as indicated by the grey horizontal 
line. 

plot the population difference about the trap origin 
AN/N = (N L - N R ) /N, where N L is the atom num- 
ber for z < 0, Ar is the atom number for z > and 
N is the initial atom number in each BSW. 

Although the BSWs start with equal populations 
in each case, an asymmetry develops over time with 
one BSW becoming increasingly populated at the 
expense of the other. This is due to a population 
transfer during each collision, as observed in [15j for 
a single BSW collision. Here the BSWs collide with 
approximate speed v — \uj r zo « 0.7mm s _1 at the 
trap centre. At this relatively large speed, the popu- 
lation transfer is small and, from [15] . we can expect 
it to be of the order of 1% N. However, the multiplic- 
ity of the collisions leads to the growth of a notice- 
able population transfer, until eventually the system 
becomes unstable against collapse. We note that the 
time at which the system becomes unstable depends 
upon the initial phase. During each collision, the rel- 
ative phase is also modified and after many collisions 
this can lead to a reversal of the population transfer. 
For example, for Atfi = 37r/4, the population differ- 
ence oscillates slowly before the system ultimately 
collapses. Since no large asymmetries are observed 
in the experimental density profiles, this validates 
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our original assumption that the relative phase must 
of either zero or 7r (or very close). Furthermore, the 
fact that all cases except it become unstable in well 
under 3 s gives further evidence that a 7r-phase dif- 
ference is crucial to maintain the dynamics of the 
two BSWs observed in the experiment. 

5. Conclusions 

In summary, we have performed a detailed anal- 
ysis of the multiple bright solitary waves (BSW) 
observed experimentally [4]. We confirm that such 
multiple BSW states are only stable if there is a 
ir phase difference between each BSW. This allows 
each BSW to contain approximately the critical 
number of atoms. Remarkably the experimental 
data implies that the atom number in each BSW 
lies consistently just under the threshold for col- 
lapse. We find that two BSWs featuring a 7r-phase 
difference undergo stable dynamics in a trap over 
long times (of order a few seconds or over 40 col- 
lisions), in excellent quantitative agreement with 
experimental measurements. In contrast, for 0- 
phase difference, the system undergoes collapse. 
For intermediate relative phases (0 < Acf> < tt) we 
observe significant population transfer between the 
BSWs and a lifetime against collapse that depends 
upon the value of the relative phase. We propose 
that these predictions could be verified by using a 
phase imprinting technique similar to that used in 
the creation of dark solitons [5] to impose a con- 
trolled relative phase on the BSWs shortly after 
their creation. 
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